grstyle init
grstyle set plain
grstyle set legend 6, nobox
grstyle set color viridis, n(8) range(0.2,0.8)

set seed 119922

* set working directory: TO BE ADJUSTED
cd ""

* load dataset
use "dataset_attr.dta", clear


******************************************************
** Attrition Analysis for Supplementary Information **
******************************************************

use "dataset_attr.dta", clear
svyset _n [pweight=inv_prob_allland]


quietly tabout attrited treat using "SI/Tables/attrition_treat.tex", replace c(freq col) style(tex) font(italic) twidth(8) f(0) clab(No. Col_%) stats(chi2)

eststo reg_att_1: quietly reg attrited base_hhead_age base_hhead_educ base_ihs_land_area base_forest_area base_any_tree_cut base_ihs_tim_price base_rent base_envprog_bi treat [pweight=inv_prob_allland], vce(cluster village_code)

eststo reg_att_2: quietly reg attrited base_hhead_age base_hhead_educ base_ihs_land_area base_forest_area base_any_tree_cut base_ihs_tim_price base_rent base_envprog_bi treat c.base_hhead_age#c.treat c.base_hhead_educ#c.treat c.base_ihs_land_area#c.treat c.base_forest_area#c.treat c.base_any_tree_cut#c.treat c.base_ihs_tim_price#c.treat c.base_rent#c.treat c.base_envprog_bi#c.treat [pweight=inv_prob_allland], vce(cluster village_code)

quietly test (c.base_hhead_age#c.treat=0) (c.base_hhead_educ#c.treat=0) (c.base_ihs_land_area#c.treat=0) (c.base_forest_area#c.treat=0) (c.base_any_tree_cut#c.treat=0) (c.base_ihs_tim_price#c.treat=0) (c.base_rent#c.treat=0) (c.base_envprog_bi#c.treat=0) 
estadd scalar test_interact = r(p): reg_att_2

* replace column allignment first column with p{6cm}

quietly esttab reg_att_1 reg_att_2 using "SI/Tables/reg_attr2.tex", l b(3) ci(3) nomtitles eqlabels(none) title("Attrition Regressions") nobase stats(test_interact N F p r2, labels("F-Test Interactions: p" "N" "F-Test" "p-value" "R2") f(%4.3f %4.0f %4.3f)) star(* 0.10 ** 0.05 *** 0.01) long replace



************************************************************************************
** Manski Bounds (imputing maximal values to 2 villages where access was refused) **
************************************************************************************

use "dataset_bounds.dta", clear
* set sampling weights to 1 for not covered village
replace inv_prob_allland = 1 if village_code == 1909
svyset _n [pweight=inv_prob_allland]

* impute outcomes opposite to results
replace soc_type_spiteful = 1 if imputed == 1
replace soc_type_generous = 1 if imputed == 1
replace soc_type_egalitarian = 0 if imputed == 1

replace for_restrictions = 1 if imputed == 1
replace for_restrictions_nfo = 1 if imputed == 1

replace index_socialnorm = 1 if imputed == 1
replace index_personalnorm = 0 if imputed == 1

* Rerun regressions from main paper

global controls i.base_subcounty strat_base_num_PFOs_vil strat_base_inc_pc_vil strat_base_dist_road_vil strat_base_l1_area_vil

quietly probit soc_type_spiteful i.treat $controls [pweight=inv_prob_allland], vce(cluster village_code)
scalar Nclust = e(N_clust)
scalar r2p = e(r2_p)
qui testparm i.treat $controls
scalar p_value = r(p)
scalar chi = r(chi2)
eststo probit_soc_spiteful_pfo: quietly margins, dydx(*) atmeans post
estadd scalar p=p_value
estadd scalar N_clust=Nclust
estadd scalar r2_p=r2p
estadd scalar chi2=chi

quietly probit soc_type_generous i.treat $controls  [pweight=inv_prob_allland], vce(cluster village_code)
scalar Nclust = e(N_clust)
scalar r2p = e(r2_p)
qui testparm i.treat $controls
scalar p_value = r(p)
scalar chi = r(chi2)
eststo probit_soc_generous_pfo: quietly margins, dydx(*) atmeans post
estadd scalar p=p_value
estadd scalar N_clust=Nclust
estadd scalar r2_p=r2p
estadd scalar chi2=chi

quietly probit soc_type_egalitarian i.treat $controls  [pweight=inv_prob_allland], vce(cluster village_code)
scalar Nclust = e(N_clust)
scalar r2p = e(r2_p)
qui testparm i.treat $controls
scalar p_value = r(p)
scalar chi = r(chi2)
eststo probit_soc_egali_pfo: quietly margins, dydx(*) atmeans post
estadd scalar p=p_value
estadd scalar N_clust=Nclust
estadd scalar r2_p=r2p
estadd scalar chi2=chi


** Forest use restrictions

eststo tobit_for_rest_for_own: quietly tobit for_restrictions i.treat $controls  [pweight=inv_prob_allland], vce(cluster village_code) ll(0) ul(1)
qui testparm i.treat $controls
estadd scalar p_value = r(p)
estadd scalar Falt = r(F)
eststo tobit_for_rest_nfo: quietly tobit for_restrictions_nfo i.treat $controls   [pweight=inv_prob_allland], vce(cluster village_code) ll(0) ul(1)
qui testparm i.treat $controls
estadd scalar p_value = r(p)
estadd scalar Falt = r(F)


** Norms

eststo tobit_pnorm_pfo: quietly tobit index_personalnorm i.treat $controls  [pweight=inv_prob_allland], vce(cluster village_code) ll(0) ul(1)
qui testparm i.treat $controls
estadd scalar p_value = r(p)
estadd scalar Falt = r(F)

eststo tobit_snorm_pfo: quietly tobit index_socialnorm i.treat $controls  [pweight=inv_prob_allland], vce(cluster village_code) ll(0) ul(1)
qui testparm i.treat $controls
estadd scalar p_value = r(p)
estadd scalar Falt = r(F)

* create table for SI
quietly esttab tobit_for_rest_for_own tobit_for_rest_nfo using "SI/Tables/reg_bounds1.tex", l b(3) ci(3) mtitles("Forest Owners" "Non-Forest Owners") title("Treatment Effect Bounds - Forest use restrictions index") eqlabels(none) drop(*.base_subcounty) nobase stats(N N_clust r2_p Falt p_value, labels("N" "Clusters" "Pseudo R2" "F-Statistic" "p") f(%4.0f %4.0f %4.3f)) star(* 0.10 ** 0.05 *** 0.01) addnotes("Stratification variables abbreviated with (strat.). All models include subcounty fixed effects." "Standard errors are clustered at the village level. Tobit models bounded between 0 and 1.") replace 

quietly esttab tobit_pnorm_pfo tobit_snorm_pfo using "SI/Tables/reg_bounds2.tex", l b(3) ci(3)  mtitles("Personal" "Social") title("Treatment Effect Bounds - Personal and Social Norms") eqlabels(none) drop(*.base_subcounty)  nobase stats(N N_clust r2_p Falt p_value, labels("N" "Clusters" "Pseudo R2" "F-Statistic" "p") f(%4.0f %4.0f %4.3f)) star(* 0.10 ** 0.05 *** 0.01) addnotes("Stratification variables abbreviated with (strat.). All models include subcounty fixed effects." "Standard errors are clustered at the village level. Tobit models bounded between 0 and 1.") replace 

quietly esttab probit_soc_spiteful_pfo probit_soc_generous_pfo probit_soc_egali_pfo using "SI/Tables/reg_bounds3.tex", l b(3) ci(3) mtitles("Spiteful" "Generous" "Egalitarian" ) title("Treatment Effect Bounds - Social Preferences") eqlabels(none) drop(*.base_subcounty)  nobase stats(N N_clust r2_p chi2 p, labels("N" "Clusters" "Pseudo R2" "Chi2 Statistic" "p") f(%4.0f %4.0f %4.3f)) star(* 0.10 ** 0.05 *** 0.01) addnotes("Stratification variables abbreviated with (strat.). All models include subcounty fixed effects." "Standard errors are clustered at the village level. Probit models reported as marginal effects.") replace 
